Preliminary SEM Analyses 2
Delta Food Web Team
10/25/2021
library(dplyr)## Warning: package 'dplyr' was built under R version 4.1.1library(psych) #for pairs.panels, but could use other packages, e.g. GGalley## Warning: package 'psych' was built under R version 4.1.1library(lavaan)## Warning: package 'lavaan' was built under R version 4.1.1library(semPlot)## Warning: package 'semPlot' was built under R version 4.1.1library(DiagrammeR)## Warning: package 'DiagrammeR' was built under R version 4.1.1Import data
combined=read.csv("data/annual_averages/annual_data_compiled.csv")
cnames=read.csv("analysis/column_names.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
source("analysis/myLavaanPlot.r")Data prep
Log transform, scale
#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
x2=x[which(!is.na(x))]
if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
else {log(x)}
}
focaldatalog = focaldata %>%
mutate_at(logvars,logtrans)
#scale data
fd=focaldatalog
fd=fd %>%
mutate_at(2:length(fvars),list("1"=lag)) %>% #lag 1
mutate_at(2:length(fvars),list("fd"=function(x) c(NA,diff(x)))) %>% #first difference
mutate_at(2:length(fvars),list("dtr"=function(x) { #detrend
x2=x
x2[x2==0]=NA
res=residuals(lm(x2~fd$year))
out=x
out[which(!is.na(x2))]=res
return(out)
})) %>%
mutate_at(-1,scale)Time series plots
Original units
Log scaled
First difference
Detrended
Bivariate plots
SEM model
Best model so far
model4='zoop=~hcope+clad+mysid
fish=~estfish+estfish_bsmt+estfish_bsot
zoop~sside+chla
chla~clams+flow
fish~zoop+flow+sside
hcope~~mysid
clad~~mysid
hcope~~clad'
modfit4=sem(model4, data=fd)
summary(modfit4, standardized=T, rsq=T)## lavaan 0.6-9 ended normally after 39 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 23
##
## Used Total
## Number of observations 40 46
##
## Model Test User Model:
##
## Test statistic 34.739
## Degrees of freedom 26
## P-value (Chi-square) 0.117
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## zoop =~
## hcope 1.000 0.607 0.615
## clad 1.341 0.334 4.013 0.000 0.814 0.913
## mysid 1.421 0.318 4.461 0.000 0.862 0.921
## fish =~
## estfish 1.000 0.758 0.766
## estfish_bsmt 1.153 0.197 5.852 0.000 0.874 0.899
## estfish_bsot 0.968 0.193 5.011 0.000 0.734 0.772
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## zoop ~
## sside -0.513 0.155 -3.298 0.001 -0.844 -0.515
## chla 0.409 0.123 3.334 0.001 0.674 0.531
## chla ~
## clams -0.415 0.112 -3.689 0.000 -0.415 -0.495
## flow -0.342 0.106 -3.215 0.001 -0.342 -0.431
## fish ~
## zoop 0.558 0.315 1.770 0.077 0.447 0.447
## flow 0.303 0.100 3.041 0.002 0.400 0.397
## sside -0.391 0.240 -1.628 0.104 -0.516 -0.315
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hcope ~~
## .mysid 0.103 0.145 0.708 0.479 0.103 0.362
## .clad ~~
## .mysid -0.097 0.154 -0.628 0.530 -0.097 -0.729
## .hcope ~~
## .clad 0.019 0.134 0.141 0.888 0.019 0.067
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hcope 0.605 0.183 3.298 0.001 0.605 0.621
## .clad 0.132 0.157 0.843 0.399 0.132 0.167
## .mysid 0.134 0.174 0.766 0.444 0.134 0.152
## .estfish 0.405 0.108 3.759 0.000 0.405 0.413
## .estfish_bsmt 0.181 0.077 2.345 0.019 0.181 0.192
## .estfish_bsot 0.365 0.098 3.729 0.000 0.365 0.404
## .chla 0.418 0.093 4.472 0.000 0.418 0.675
## .zoop 0.123 0.103 1.197 0.231 0.334 0.334
## .fish 0.171 0.076 2.235 0.025 0.297 0.297
##
## R-Square:
## Estimate
## hcope 0.379
## clad 0.833
## mysid 0.848
## estfish 0.587
## estfish_bsmt 0.808
## estfish_bsot 0.596
## chla 0.325
## zoop 0.666
## fish 0.703semPaths(modfit4, "std", edge.label.cex = 1, residuals = F)
semPaths(modfit4, "par", edge.label.cex = 1, residuals = F)
# lavaan::parameterestimates(modfit4)
labels1 <- createLabels(modfit4, cnames)
## Diagram without model covariances:
myLavaanPlot(model=modfit4, labels=labels1,
node_options=list(shape="box", fontname="Helvetica"),
coefs=TRUE, stand=TRUE, covs=FALSE, sig=0.05,
width=c("regress","latent"),
color=c("regress","latent"))## Diagram with model covariances:
myLavaanPlot(model=modfit4, labels=labels1,
node_options=list(shape="box", fontname="Helvetica"),
coefs=TRUE, stand=TRUE, covs=TRUE, sig=0.05,
width=c("regress","latent","covs"),
color=c("regress","latent","covs"))Same model with detrended data
model4b='zoop=~hcope_dtr+clad_dtr+mysid_dtr
fish=~estfish_dtr+estfish_bsmt_dtr+estfish_bsot_dtr
zoop~sside_dtr+chla_dtr
chla_dtr~clams_dtr+flow_dtr
fish~zoop+flow_dtr+sside_dtr'
modfit4b=sem(model4b, data=fd)
summary(modfit4b, standardized=T, rsq=T)## lavaan 0.6-9 ended normally after 29 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 20
##
## Used Total
## Number of observations 40 46
##
## Model Test User Model:
##
## Test statistic 34.591
## Degrees of freedom 29
## P-value (Chi-square) 0.218
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## zoop =~
## hcope_dtr 1.000 0.743 0.754
## clad_dtr 0.867 0.245 3.543 0.000 0.644 0.641
## mysid_dtr 1.072 0.269 3.990 0.000 0.796 0.795
## fish =~
## estfish_dtr 1.000 0.441 0.485
## estfsh_bsmt_dt 2.169 0.769 2.818 0.005 0.956 0.971
## estfsh_bst_dtr 1.296 0.482 2.688 0.007 0.571 0.580
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## zoop ~
## sside_dtr -0.303 0.182 -1.670 0.095 -0.408 -0.272
## chla_dtr 0.390 0.138 2.832 0.005 0.525 0.492
## chla_dtr ~
## clams_dtr 0.056 0.128 0.435 0.664 0.056 0.063
## flow_dtr -0.411 0.137 -3.007 0.003 -0.411 -0.435
## fish ~
## zoop 0.104 0.100 1.036 0.300 0.175 0.175
## flow_dtr 0.248 0.106 2.348 0.019 0.563 0.558
## sside_dtr -0.093 0.099 -0.946 0.344 -0.212 -0.141
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hcope_dtr 0.419 0.145 2.892 0.004 0.419 0.431
## .clad_dtr 0.595 0.161 3.685 0.000 0.595 0.590
## .mysid_dtr 0.369 0.149 2.467 0.014 0.369 0.368
## .estfish_dtr 0.633 0.148 4.267 0.000 0.633 0.765
## .estfsh_bsmt_dt 0.055 0.191 0.291 0.771 0.055 0.057
## .estfsh_bst_dtr 0.644 0.160 4.017 0.000 0.644 0.664
## .chla_dtr 0.697 0.156 4.472 0.000 0.697 0.795
## .zoop 0.395 0.168 2.352 0.019 0.716 0.716
## .fish 0.121 0.077 1.573 0.116 0.620 0.620
##
## R-Square:
## Estimate
## hcope_dtr 0.569
## clad_dtr 0.410
## mysid_dtr 0.632
## estfish_dtr 0.235
## estfsh_bsmt_dt 0.943
## estfsh_bst_dtr 0.336
## chla_dtr 0.205
## zoop 0.284
## fish 0.380semPaths(modfit4b, "std", edge.label.cex = 1, residuals = F)
semPaths(modfit4b, "par", edge.label.cex = 1, residuals = F)
Old models (not run)
model1='clams~year+flow
sside~year+flow
chla~flow+clams #+year?
hcope~chla+pcope #+year?
clad~chla+pcope
#pcope~year+clams
mysid~hcope+pcope
estfish~hcope+pcope+mysid+sside'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)
residuals(modfit1)
model1='clams~year+flow
sside~year+flow
foodsupply=~chla+hcope+clad+mysid
foodsupply~clams+flow+sside
fish=~estfish+estfish_bsmt
fish~foodsupply'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)
model1='foodsupply=~chla+hcope+clad+mysid
fish=~estfish+estfish_bsmt
foodsupply~clams+flow+sside
fish~foodsupply'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)
model2='foodsupply=~chla+hcope+clad+mysid
fish=~estfish+estfish_bsmt
foodsupply~clams+flow+sside
fish~foodsupply+flow+sside'
modfit2=sem(model2, data=fd)
summary(modfit2, standardized=T, rsq=T)
semPaths(modfit2, "std", edge.label.cex = 1, residuals = F)
residuals(modfit2)
model3='zoop=~hcope+clad+mysid
fish=~estfish+estfish_bsmt
zoop~clams+flow+sside+chla
chla~clams #+tphos+nitrate+ammonia
fish~zoop+flow+sside'
modfit3=sem(model3, data=fd)
summary(modfit3, standardized=T, rsq=T)
semPaths(modfit3, "std", edge.label.cex = 1, residuals = F)
residuals(modfit3, type="cor")
modificationindices(modfit3)
model4='zoop=~hcope+clad+mysid
fish=~estfish+estfish_bsmt+estfish_bsot
zoop~sside+chla
chla~clams+flow
fish~zoop+flow+sside
hcope~~mysid
clad~~mysid
hcope~~clad'
modfit4=sem(model4, data=fd)
summary(modfit4, standardized=T, rsq=T)
semPaths(modfit4, "std", edge.label.cex = 1, residuals = F)
semPaths(modfit4, "par", edge.label.cex = 1, residuals = F)
residuals(modfit4, type="cor")
modificationindices(modfit4)#using estfish instead of smelt
model1f='chla~year+flow+secchi
hcope~year+flow+temp+chla+pcope+secchi
pcope~year
mysid~year+hcope+pcope
estfish~year+secchi+hcope+pcope+mysid'
modfit1f=sem(model1f, data=fd)
#summary(modfit1f, standardized=T, rsq=T)
semPaths(modfit1f, "std", edge.label.cex = 1, residuals = F)
#using marfish instead of smelt
model1g='chla~year+flow+secchi
hcope~year+flow+temp+chla+pcope+secchi
pcope~year
mysid~year+hcope+pcope
marfish~year+secchi+hcope+pcope+mysid'
modfit1g=sem(model1g, data=fd)
#summary(modfit1g, standardized=T, rsq=T)
semPaths(modfit1g, "std", edge.label.cex = 1, residuals = F)Detrended (anomalies) without year effect
#smelt
model2a='chla_dtr~flow_dtr+temp_dtr+secchi_dtr
hzoop_dtr~flow_dtr+temp_dtr+secchi_dtr+chla_dtr+pcope_dtr
pcope_dtr~flow_dtr+temp_dtr+secchi_dtr
mysid_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr
smelt_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr+mysid_dtr'
modfit2a=sem(model2a, data=fd)
#summary(modfit2a, standardized=T, rsq=T)
semPaths(modfit2a, "std", edge.label.cex = 1, residuals = F)
#estuarine fishes
model2b='chla_dtr~flow_dtr+temp_dtr+secchi_dtr
hzoop_dtr~flow_dtr+temp_dtr+secchi_dtr+chla_dtr+pcope_dtr
pcope_dtr~flow_dtr+temp_dtr+secchi_dtr
mysid_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr
estfish_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr+mysid_dtr'
modfit2b=sem(model2b, data=fd)
#summary(modfit2b, standardized=T, rsq=T)
semPaths(modfit2b, "std", edge.label.cex = 1, residuals = F)